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C^ Abstract 



We propose a solution to the image deconvolution problem where the convolution kernel or point spread function 

(PSF) is assumed to be only partially known. Small perturbations generated from the model are exploited to produce a 

few principal components explaining the PSF uncertainty in a high dimensional space. Unlike recent developments on 

(~| blind deconvolution of natural images, we assume the image is sparse in the pixel basis, a natural sparsity arising in 

I magnetic resonance force microscopy (MRFM). Our approach adopts a Bayesian Metropolis-within-Gibbs sampling 



I. Introduction 



framework. The performance of our Bayesian semi-blind algorithm for sparse images is superior to previously proposed 
semi-blind algorithms such as the alternating minimization (AM) algorithm and blind algorithms developed for natural 
jy! images. We illustrate our myopic algorithm on real MRFM tobacco virus data. 
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C^ Recently, a new 3D imaging technology called magnetic resonance force microscopy (MRFM) has been developed. 

pg The principles of MRFM were introduced by Sidles IT], ||2l, ^ who described its potential for achieving 3D atomic 

. . scale resolution. In 1992 and 1996, Rugar et al. p], fS) reported experiments that demonstrated the feasibility of 

> 

MRFM and produced the first MRFM images. More recently, MRFM volumetric spatial resolutions of less than 

lOnm have been demonstrated for imaging a biological sample 0. The signal provided by MRFM is a so- 
called force map that is the 3D convolution of the atomic spin distribution and the point spread function (PSF) 
Q. This formulation casts the estimation of the spin density from the force map as an inverse problem. Several 
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approaches have been proposed to solve this inverse problem, i.e., to reconstruct the unknown image from the 
measured force map. Basic algorithms rely on Wiener filters (SI, ||9l, whereas others are based on iterative 
least squares reconstruction approaches Q, lID, ifTOl . More recently, this problem has been addressed within the 
Bayesian estimation framework ifTTI . lfT2l . 

However, all of these reconstruction techniques require prior knowledge of the device response, namely the PSF. 
As shown by Mamin et al. lfT3l . this PSF is a function of several parameters specified by the physics of the device 
including: mass of cantilever probe, ferromagnetic constant of probe tip, and external field strength. Unfortunately, 
in practice the physical parameters that tune the response of the MRFM tip are only partially known. In such 
circumstances, the PSF used in the reconstruction algorithm is mismatched to the true PSF of the microscope 
and the quality of standard MRFM image reconstruction will suffer if one does not account for this mismatch. 
Estimating the unknown image and the partially known PSF jointly is usually referred to as semi-blind lfT4l . ifTSJI 
or myopic lfT6l . ifTTJI deconvolution, and this is the approach taken in this paper. The myopic image restoration 
problem was previously studied within a hierarchical Bayesian framework lITSl with partially known blur functions 
in many applications, including natural and astronomical imaging lfT9l . Il20l . This previous work lfT9l . Il20l models 
the deviation of the PSF as uncorrected zero mean Gaussian noise. The authors of ||2T| considered an extension 
of this model to a non-sparse, simultaneous autoregression (SAR) prior model for both the image and point spread 
function. Compared to ^T\\, recent papers on single motion deblurring in photography ll22l . ||23]| use heavier tailed 
distributions for the gradient of images and an exponential distribution for the PSF. In addition, the algorithm in |l22l 
separately identifies the PSF using a multi-scale approach to perform conventional image restoration. The authors 
of II23I proposed an image prior to reduce ringing artifacts from blind deconvolution of photo images. This paper 
considers an alternative model, appropriate to MRFM but significantly different from photography, that imposes 
sparsity on the image and an empirical Bayes prior on the PSF. 

To mitigate the effects of PSF mismatch on MRFM sparse image reconstruction, a non-Bayesian alternating 
minimization (AM) algorithm 1241 was proposed by Herrity et al. which showed robust performance. In this paper, 
we propose a hierarchical Bayesian approach to semi-blind image deconvolution that exploits prior information on 
the PSF model. This is a semi-blind modification of the Bayesian MRFM reconstruction approach of Dobigeon et al. 
lfT2l that uses an adaptive tuning scheme to produce a Bayesian estimate of the PSF and a Bayesian reconstruction 
of the image. The contribution of this paper is a novel Bayesian approach to a joint estimation of PSFs and images. 
We represent the PSF on a truncated orthogonal basis, where the basis elements are the singular vectors in the 
singular value decomposition of the family of perturbed nominal PSFs. A Gaussian prior model specifies a log 
quadratic Bayes prior on deviations from the nominal PSF. Our approach is related to the recent papers of Tzikas 
et al. f25l and Orieux et al. [26]. In [25 1 a pixel-wise, space-invariant Gaussian kernel basis is assumed with a 
gradient based image prior Orieux et al. introduced a Metropolis -within-Gibbs algorithm to estimate the parameters 
that tune the device response. The strategy II26I focuses on reconstruction with smoothness constraints and requires 
recomputation of the entire PSF at each step of the algorithm. This is computationally expensive, especially for 
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complex PSF models such as in the MRFM instrument. Here, we propose an alternative that consists of estimating 
the deviation from a given nominal PSF. More precisely, the nominal point response of the device is assumed 
known and the true PSF is modeled as a small perturbation about the nominal response. Since we only need to 
estimate linear perturbations about the nominal PSF relative to a low dimensional precomputed and truncated basis 
set, this leads to reduction in computational complexity and an improvement in convergence as compared to ESl 
and ll26l . We approximate the full posterior distribution of the PSF and the image using samples generated by a 
Markov Chain Monte Carlo algorithm. Simulations and comparisons to other state-of-the-art blind deconvolution 
algorithms are presented and quantify the advantages of our algorithm for myopic sparse image reconstruction. We 
then apply it to real MRPM tobacco virus data made available by our IBM collaborators. 

This paper is organized as follows: Section [H] formulates the problem. Section [III] covers the Bayesian framework 
of image modeling and the following Section IV proposes a solution in this framework. Section [V] shows simulation 



results and an application to the real MRFM data. 



II. Forward Imaging and PSF Model 

Let X denote the /i x . . . x ?„ unknown ri-D positive spin density image to be recovered (e.g., n = 2 or 71 = 3) 
and X e R*^ denote the vectorized version of X with AI — I1I2 . ■ -In- This image is to be reconstructed from a 
collection of P(= M) measurements y — [yi, . . . , j/p] via the following noisy transformation; 

y = rKx) + n, (1) 

where T {■, •) is the 71-dimensional convolution operator or the mean response function E[y|K;,x], n is a P x 1 
observation noise vector and k is the kernel modeling the response of the imaging device. 

A typical PSF for MRFM is shown in Mamin et al. f[3l for horizontal and vertical MRFM tip configurations. 
In ([T]) n is an additive Gaussian noise sequence, independent of x, distributed according to n ^ Af (O, c^Ip). The 
PSF is assumed to be known up to a perturbation Ak about a known nominal Kq: 



K — Kq + Ak. (2) 

In the MRFM application the PSF is described by an approximate parametric function that depends on the 
experimental setup. Based on the physical parameters (gathered in the vector Q tuned during the experiment 
(external magnetic field, mass of the probe, etc.), an approximation Kq of the PSF can be derived. However, due 
to model mismatch and experimental errors, the true PSF k may deviate from the nominal PSF Kq. 

If a vector of the nominal values of parameters ^q for the parametric PSF model Kgen(C) is known, then direct 
estimation of a parameter deviation, A^, can be performed by evaluation of Hgen{Co + ^C)- However, in MRFM, 
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as shown by Mamin et al. ifTSJI . HgeniC) is a nonlinear function with many parameters that are required to satisfy 
'resonance conditions' to produce a meaningful MRFM PSF. This makes direct estimation of the PSF difficult. 



In this paper, we take a similar approach to the 'basis expansions' in 1271 Chap. 5], 1251 for approximation of the 
PSF deviation A/« as linear models, this deviation is that Ak can be expressed as a linear combination of elements 
of an a priori known basis Wk, k — 1, . . . ,K, 

K 

AK = ^AfcVfc, (3) 

fc=i 
where {vfe}j,^j^ a: is a set of basis functions for the PSF perturbations and A^, k = 1,...,K are unknown 

coefficients. To emphasize the influence of these coefficients on the actual PSF, n will be denoted k (A) with 

A = [Al, . . . , Xk] ■ With these notations, ^ can be rewritten: 

y = r(/^(A),x) + n = H(A)x + n, (4) 

where H (A) is an P x M matrix that describes convolution by the PSF kernel k (A). 

We next address the problem of estimating the unobserved image x and the PSF perturbation An under sparsity 
constraints given the measurement y and the biUnear function T {■,■). 

III. Hierarchical Bayesian model 

A. Likelihood function 

Under the hypothesis that the noise in ([TJ is Gaussian, the observation model UkelUiood function takes the form 

where ||-|| denotes the standard ^2 norm: ||x|| = x-'^x. This likelihood function will be denoted f{y\9), where 
0={x,A,a2}. 

B. Parameter prior distributions 

In this section, we introduce prior distributions for the parameters 9. 

1) Image prior: As the prior distribution for Xi, we adopt a mixture of a mass at zero and a single-sided 

exponential distribution: 

f {xi\w, a) = (1 - w)S (xi) H exp -] 1r. (xi) , (6) 

a V a / + 

where we [0, 1], a e [0, oo), S (•) is the Dirac function, R^ is a set of real open interval (0, oo) and 1^ (x) is the 

indicator function of the set E: 

' 1, if X e E, 
1e (x) = <( (7) 

0, otherwise. 
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By assuming the components Xi to be a conditionally independent (i ~ l,...,Af) given w,a, the following 
conditional prior distribution is obtained for the image x: 

M 

This image prior is similar to the LAZE distribution (weighted average of a Laplacian probability density function 
(pdf) and an atom at zero) used, for example, in Ting et al. ESJI . ifTTI . As motivated by Ting et al. and Dobigeon 
et al. ifTTIl . lfT2l . the image prior in (6) has the interesting property of enforcing some pixel values to be zero, 
reflecting the natural sparsity of the MRFM images. Furthermore, the proposed prior in (|6]l ensures positivity of 
the pixel values (spin density) to be estimated. 

2} PSF parameter prior: We assume that the parameters Ai, . . . , \k are a priori independent and uniformly 
distributed over known intervals associated with error tolerances centered at 0. Specifically, define the interval 



5fe = hAAfc,AAfc] (9) 



and define the distribution of A as 



/w=n^i^.(^^)' (10) 

k=l 

with A = [Ai,...,A/f] . In our experiment, AA^'s are set to be large enough to be non-informative, i.e., an 
improper, flat prior. 

3) Noise variance prior: A non-informative Jeffreys' prior is selected as prior distribution for the noise variance: 

f{^')^^2 (11) 

This model is equivalent to an inverse gamma prior with a non-informative Jeffreys' hyperprior, which can be seen 
by integrating out the variable corresponding to the hyperprior llT2l . 

C. Hyperparameter priors 

Define the hyperparameter vector associated with the image and noise variance prior distributions as $ = {a, w}. 
In our hierarchical Bayesian framework, the estimation of these hyperparameters requires prior distributions in the 
hyperparameters. These priors are defined in lfT2l but for completeness their definitions are reproduced below. 

1) Hyperparameter a: A conjugate inverse-Gamma distribution is assumed for hyperparameter a: 

a\a r^ig {ao,ai) , (12) 

with a = [aQ,ai] . The fixed hyperparameters ao and ai have been chosen to produce a vague prior, i.e.. 



ao = ai = 10 



-10 
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2) Hyperparameter w: A uniform distribution on the simplex [0, 1] is selected as prior distribution for the mean 
proportion of non-zero pixels: 

w^U {[0,1]). (13) 

Assuming that the individual hyperparameters are independent the full hyperparameter prior distribution for $ 
can be expressed as: 

/($|a) = /H/(a) 

1 a '^^^^ 

°^ ^j^;^ e^p (-^) i[o,i] (^) 1r+ («) ' 

D. Posterior distribution 

The posterior distribution of {0, $} is: 

/(0,*|y)cx/(y|0)/(0|*)/(*), (15) 

with 

/(0|*) = /(x|a,^i;)/(A)/(a2), (16) 

where / {y\9) and / ($) have been defined in (|5]) and ( [T4] i. The conjugacy of priors in this hierarchical structure 
allows one to integrate out the parameters a^, and the hyperparameter $ in the full posterior distribution ([TSj, 
yielding: 

/ (x, A|y, ao, ai) ex / / (0, *|y) dwdada^ ex (17) 



Se (1 + rii, 1 + no) T {ni + ao) 



Y ] ^ 1 



||y - T (a. (A) , x)r (ll^lli + «i)"^+"'' tA 2AA, 
where Be is the beta function and F is the gamma function. 

The next section presents the Metropolis-within-Gibbs algorithm 1291 that generates samples distributed according 
to the posterior distribution / (x, A|y). These samples are then used to estimate x and A. 

IV. Metropolis-within-Gibbs algorithm 

FOR semi-blind SPARSE IMAGE RECONSTRUCTION 

We describe in this section a Metropolis-within-Gibbs sampling strategy that allows one to generate samples 
< x^*), a'*' > distributed according to the posterior distribution in (TT) . As sampling directly from ^TT) is a 

difficult task, we will instead generate samples distributed according to the joint posterior / (x, A, cr^|y,ao, ai). 
Sampling from this posterior distribution is done by alternatively sampling one of x, A, a^ conditioned on all other 
variables lf30l . Ifl2l . The contribution of this work to |12| is to present an algorithm that simultaneously estimates 
both the image and PSF. The algorithm results in consistently stable output images and PSFs. 
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The main steps of our proposed sampling algorithm are given in subsections IV-A through IV-C (see also 
Algorithm [l]|. 

Algorithm 1 Metropolis-within-Gibbs sampling algorithm for semi-blind sparse image reconstruction 

1: % Initialization: 

2: Sample the unknown image x*^*^' from pdf in (|8]l, 
3: Sample the noise variance a^lo) fj-om the pdf in ([TTJ, 

4: % Iterations: 

5: for t= 1,2,..., do 

6: Sample hyperparameter w'*^ from the pdf in ( [T9| l, 

7: Sample hyperparameter a'*) from the pdf in ( |20] i, 

8: For i — I, . . . , M, sample the pixel intensity xl from the pdf in pT| ), 

9: For k = 1, . . . ,K, sample the PSF parameter X^. from the pdf in ( |23] l, by using Metropolis-Hastings step 

(see Algo. |2]), 
10: Sample the noise variance ct^*^*^ from the pdf in (|26|, 

11: end for 



A. Generation of samples according to f (x A, cr^, y, ao, ai ) 

To generate samples distributed according to / (x |A, (7^,y, ao, ai), it is convenient to sample according to 
/ (x, w, a I A, cr^, y, aoi cki ) by the following 3-step procedure. 

1) Generation of samples according to f {w |x, ao, ai ).• The conditional posterior distribution of w is 

fiw\x)(x{l-wr°+^-^w"'+^-\ (18) 

where ni = |1x||q and no = M — |lx||g. Therefore, generation of samples according to /(w|x) is achieved as 
follows: 

w|x-6e(l + ni,l + no). (19) 

2) Generation of samples according to f (a |x).- The joint posterior distribution ( [T5| ) yields: 

a|x,ao,ai -- ZtJ (||x||(, + ao, ||x||^ + ai) . (20) 

3) Generation of samples according to / (x Iw, a, A^cr^, y ).• The posterior distribution of each component Xi 
(i = 1, . . . , M) given all other variables is derived as: 

/ (xj|w,a, A,cr^,x_i,y) ex (1 ~ Wi)S (xi) + Wi(l)+ {xi\iii,rjl) , (21) 
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where x_j stands for the vector x whose ith component has been removed and /i^ and rjf are defined as the 

following: 

2 ^' 2/^hfe, 1\ 

with hi and e^ defined in Appendix [A] 

In pT| ), (f>+ (•, m, s^) stands for the pdf of the truncated Gaussian distribution defined on R^ with hidden mean 
m and hidden variance s^. Therefore, from ( |2T] l, Xi\w,a,X(T'^,x_i,y is a Bernoulli truncated-Gaussian variable 
with parameter (w^ ,iii,rif). 

To summarize, generation of samples distributed according to f (x\w, a"^, a, , y ) can be performed by updating 
the coordinates of x using M Gibbs moves (requiring generation of Bernoulli truncated-Gaussian variables). A 
pixel-wise fast and recursive sampling strategy is presented in Appendix IA] and an accelerated sparsity enforcing 
simulation scheme is described in Appendix IB] 

B. Generation of samples according to / (A x,cr^,y) 

The posterior distribution of the parameter Aj. conditioned on the unknown image x, the noise variance a^ and 
the other PSF parameters {Aj} ,j, is 

||y-r(«(A),x)|f 



/(Afc|A_fc,x,cr ,y) cxcxp 



I5. (Afc) , (23) 



2a2 

with A_fc = {^j} -Lk- ^^ summarize in Algorithm l2| a procedure for generating samples distributed according 
to the posterior in ( |23] l using a Metropolis-Hastings step with a random walk proposition ||29l from a centered 
Gaussian distribution. In order to sample efficiently, the detailed procedure of how to choose an appropriate value 
of the variance sj. of the proposal distribution for Afc is described in Appendix O At iteration t of the algorithm, 
the acceptance probability of a proposed state A^ is: 

^A<*)->A* "" ™i^ (1' "fcl-Sfo (^*k)) , (24) 

with 



2a log afc = 



y-T(/.(Ai*0,x^ ' 



-\\y-T{K{Xl),^)\\\ (25) 



Computing the transformation T(-,) at each step of the sampler can be computationally costly. Appendix [a| 
provides a recursive strategy to efficiently sample according to / (A Ix, cr'^,y). 

C. Generation of samples according to f (a^ |x,y, A) 

Samples (ct^)^*^ are generated according to the inverse gamma posterior 

/,.^|x,,.A)=is(f.fcl(|at^). a« 
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Algorithm 2 Sampling according to / (Afc|A_s;,x, i7^,y) 



1: Sample e-7V(0,s2), 

2: Propose AJ. according to A^. = A). + £, 

3: Draw u^^U {[0,1]), 



4: Set A 



(*+i) _ 



XI, ifu<pit 
Aj, , otherwise. 



xAf 



where U (E) stands for the uniform distribution on the set E. 



V. Experiments 



In this section we present simulation results that compare the proposed semi-blind Bayesian deconvolution 
algorithms of Section IV with the non-blind method fT2|, the AM algorithm ll24l . and other blind deconvolution 
methods. Here a nominal PSF Kq was selected that corresponds to the mathematical MRFM point response model 
proposed by Mamin et al. ifTSJI . 



Using our MCMC algorithm described in Sec. IV the MMSE estimators of image and PSF are approximated by 
ensemble averages over the generated samples after the burn-in period. The joint MAP estimator is selected among 
the drawn samples, after the stationary distribution is achieved, such that it maximizes the posterior distribution 

/(x,A|y) 131 1. 



TABLE I 

Parameters used to compute the MRFM PSF. 



Parameter 




Value 






Description 


Name 




Amplitude of external magnetic field 


Best 


9.4 X 10^ G 


Value of -Bmag in the resonant slice 


Bres 


1.0 X 10-* G 


Radius of tip 


i?o 


4.0 nm 


Distance from tip to sample 


do 


6.0 nm 


Cantilever tip moment 


m 


4.6 X IQS emu 


Peak cantilever oscillation oscillation 


Xpii 


0.8 nm 


Maximum magnetic field gradient 


^max 


125 
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A. Simulation on simulated sparse images 



We performed simulations of MRFM measurements for PSF and image models similar to those described in 
Dobigeon et al. {V2\. The signal-to-noise ratio was set to SNR = lOdB. Several 32 x 32 synthetic sparse images, 



one of which is depicted in Fig. 1(a) were used to produce the data and were estimated using the proposed Bayesian 
method. The assumed PSF /to is generated following the physical model described in Mamin et al. fTS]] when the 
physical parameters are tuned to the values displayed in Table |l] This yields a 11 x 11 2-dimensional convolution 



kernel represented in Fig. 2(a) We assume that the true PSF n comes from the same physical model where the 
radius of the tip and the distance from tip to sample have been mis-specified with 2% error as R = Rq — 2% = 3.92 



and d = da + 2% = 6.12. This leads to the convolution kernel depicted in Fig. 2(b) The observed measurements 



y, shown Fig. 1(b) are a 32 x 32 image of size P = 1024. 




5 20 25 30 



^0 




(a) Sparse true image (||x||o 

11) 



(b) Raw MRFM observation 



Fig. 1 . Simulated true image and MRFM raw image exhibiting superposition of point responses (see Fig. [2l and noise. 





(a) Nominal MRFM PSF 
Fig. 2. Assumed PSF kq and actual PSF k. 



(b) True MRFM PSF 



The proposed algorithm requires the definition of K basis vectors v^, k = l,...,K, that span a subspace 
representing possible perturbations An. We empirically determined this basis using the following PSF variational 
eigendecomposition approach. A set of 5000 experimental PSFs iij, j ~ 1, ■ • • , 5000, were generated following the 
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model described in Mamin et al. ifTSl with parameters d and R randomly drawn according to Gaussian distributionH 
centered at the nominal values cLq, Rq, respectively. Then a standard principal component analysis (PCA) of the 
residuals {kj — Kq} .^^ ^^^^^ , by allowing the maximum variance over the parameters that produce valid MRFM 
PSFs, was used to identify K = 4 principal axes that are associated with the basis vectors vj.. The necessary 
number of basis vectors, K ^ 4 here, was determined empirically by detecting a knee at the scree plot shown 
in Fig. [3] The first four eigenfunctions, corresponding to the first four largest eigenvalues, explain 98.69% of the 
observed perturbations. The 4 principal patterns of basis vectors are depicted in Fig. |4] 




n eigenvalues 



Fig. 3. Scree plot of residual PCA approximation error in Z2 norm (magnitude is normalized up to the largest value, i.e. A„ 



The proposed semi-blind Bayesian reconstruction algorithm was applied to estimate both the sparse image and 
the PSF coefficients of v^'s, using the prior in (|6]l. From the observation shown in Fig. 1(b) the PSF estimated 
by the proposed algorithm is shown in Fig. |5(a)| and is in good agreement with the true one. The corresponding 



maximum a posteriori estimate (MAP) of the unknown image is depicted in Fig. 6(a) The obtained coefficients of 
the PSF-eigenfunctions are close to true coefficients (Fig. |7|. 



' We used a PSF generator provided by Dan Rugar's group at IBM 1 13|. The variances of the Gaussian distributions are carefully tuned so 
that their standard deviations produce a minimal volume ellipsoid that contains the set of valid PSF's of the form specified in 1 13 1. 



March 1, 2013 



DRAFT 



12 





2 4 B B 10 



2 4 6 B 10 
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2 4 B B 10 




2 4 6 



Fig. 4. The K = 4 principal vectors (vfe) of the perturbed PSF, identified by PCA. 



B. Comparison to Other Methods 



For comparison to a non-blind method. Fig. 6(b) shows the estimate using the Bayesian non-bhnd method lfT2ll 



with a mismatched PSF. Fig. 6(c) shows the estimate generated by the AM algorithm. The nominal PSF described 



in Section V-A is used in the AM algorithm and hereafter in other semi-blind algorithms, and the parameter values 
of AM algorithm were set empirically according to the procedure in [24J. Our proposed algorithm appears to 
outperform the others (Fig. l6]l while preserving fast convergence (Fig. iTJi. 

Quantitative comparisons were obtained by generating different noises in 100 independent trials for a fixed true 
image. Here, six different true images with six corresponding different sparsity levels (|lx||o = 6, 11, 18, 30, 59, 97) 
were tested. Fig. [8] presents the two histograms of the results with the six sets in the corresponding two error 
criteria, j|x — x|p, ||x||o, respectively, both of which indicate that our method performs better and is more stable 
than the other two methods. 

Fig. l9] shows reconstruction error performance for several measures of error used in Ting et al. fTTl and Dobigeon 
et al. 1 12 1 to compare different reconstruction algorithms for sparse MRFM images. Notably, compared to the AM 
algorithm that aims to compensate 'blindness' of the unknown PSF and the previous non-blind method, our method 
reveals a significant performance gain under most of the investigated performance criteria and sparsity conditions. 
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(a) Proposed method 



(b) Amizic's method 



(c) Almeida's method 



(d) Tzikas' method 



Fig. 5. Estimated PSF k of MRFM tip using our semi-bUnd method, Amizic's method (using TV prior), Almeida's method (using progressive 
regularization), and Tzikas' method (using the kernel basis PSF model), respectively. For fairness, we used sparse image priors for the methods. 
(See Section [V-B I for details on the methods.) 





20 25 30 




^0 




5 10 15 20 25 30 



(a) MAP, proposed method 



(b) MAP, Bayesian non-blind 
method with kq 



(c) AM 



(d) Almeida's method 




10 15 20 25 30 



(e) Tzikas' method 
Fig. 6. The true sparse image and estimated images from Bayesian non-blind, AM, our semi-blind, Almeida's, and Tzikas' methods. 



In addition to the AM and non-blind comparisons shown in Fig. [81 we made direct comparisons between our 
sparse MRFM reconstruction method and several state-of-the-art blind image reconstruction methods li32l . Il33l . 
Il25l . II22I . II23I . In all cases the algorithms were initialized with the nominal, mismatched PSF and were applied 
to a sparse MRFM- type image Uke in Fig. [T] For a fair comparison, we made a sparse prior modification in the 
image model of other algorithms. The total variation (TV) based prior for the PSF suggested by Amizic et al. Il32l 
was also implemented. The obtained PSF from this method was considerably worse than the one estimated by our 
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Fig. 7. Estimated PSF coefficients for 4 PCs over 200 iterations. These curves sliow fast convergence of our algoritlim. 'Ideal coefficients' are 
the projection values of the true PSF onto the space spanned by four principal PSF basis. 
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(a) Histograms of the normalized I2 norm en'ors. x-axis is | 
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Fig. 8. Histograms of I2 and Zq norm of the reconstruction error. Note that the proposed semi-blind reconstructions exhibit smaller mean error 
and more concentrated error distribution than the non-blind method of 1 12iJ and the alternating minimization method of 1241 . 



proposed method (see Fig. 5(b) 1 resulting in a very poor quality image deconvolutiorp] 



The recent blind deconvolution method proposed by Almeida et al. |33| utilizes the 'sharp edge' property in 
natural images, with initial, high regularization in order to effectively evaluate the PSF. This iterative approach 
has a sequentially decreasing regularization parameter to reconstruct fine details of the image. Adapted to sparse 
images, this method performs worse than our method, in terms of image and PSF estimation errors. The PSF and 



image estimates from Almeida's method are presented in Fig. 5(c) and 6(d) respectively. 



- Because this PSF is wrongly estimated and similar to the 2D delta function, the image estimate looks similar to the denoised version of 
observation, so we omit the image estimate. 
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Fig. 9. Error bar graphs of results from our myopic deconvolution algorithm. For several image x's of different sparsity levels, errors are 
illustrated with standard deviations. (Some of the sparsity measure and residual errors are too large to be plotted together with results from 
other algorithms.) 



Tzikas et al. f25l uses a similar PSF model to our method using basis kernels. However, no sparse image prior 
was assumed in |25 1 making it unsuitable for sparse reconstruction problems such as the MRFM problem considered 
in the paper. For a fair comparison, we applied the suggested PSF model lIZSJI along with our sparse image prior 
The results of PSF and image estimation and the performance graph are shown in Fig. 5(d) Fig. 6(e) and Fig. l9] 
respectively. In terms of PSF estimation error, our algorithm outperforms the others. 

We also compared against the mixture model -based algorithm of Fergus et al. Il22l . and the related method of 
Shan et al. |23|, which are proposed for blind deconvolution of shaking/motion blurs and do not incorporate any 
sparsity penalization. When applied to the sparse MRFM reconstruction problem the algorithms of |22| and ||23]| 
performed very poorly (produced divergent or trivial solutions; not shown due to space limitations). This poor 
performance is likely due to the fact that the shape of the MRFM PSF and the sparse image model are significantly 
different from those in blind deconvolution of camera shaking/motion blurs. The generalized PSF model of ll22l 
and II23I with the sparse image prior is Tzikas' model ||25]| . which is described above. 



We used the Iterative Shrinkage/Thresholding (1ST) algorithm with a true PSF to lower bound our myopic 
reconstruction algorithm. The 1ST algorithm effectively reconstructs images with a sparsity constraint ||34| . From 
Fig. 9(b) the performance of our result is as good as that of the oracle 1ST. In Table [ll] we present comparisons of 
computation timajof the proposed sparse semi-blind Bayes reconstruction to that of several other algorithms. 



'Matlab is used under Windows 7 Enterprise and HP-Z200 (Quad 2.66 GHz) platform. 
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TABLE II 

Computation time of algorithms (in seconds), for the data in Fig.[T] 



Proposed method 


19.06 


Bayesian nonblind fTl] 


3.61 


1ST (34) 


0.09 


AM |24| 


0.40 


Almeida's method (33| 


5.63 


Amizic's method 1321 


5.69 


Tzikas' method 1251 


20.31 



C. Application to 3D MRFM image reconstruction 



In this section, we apply the semi-blind Bayesian reconstruction algorithm to the 3D MRFM tobacco virus data 
||6l shown in Fig. 10 The necessary modification for our algorithm to apply to 3D data is simple because the 
extension of our 2D pixel-wise sampling method requires only one more added dimension to extend to 3D basis 
vectors and 3D convolution kernel. As seen in Appendix |AJ the voxel-wise update of a vectorized image x can be 
generalized to nD data. This scalability is another benefit of our algorithm. The implementation of AM algorithm 
is impractical due to its slow convergence rates 1241 . Here we only consider Bayesian methods. The additive noise 
is assumed Gaussian consistently with ||4|, (|6l, so the noise model in paragraph |III-A| is applied here. 





z=62 nm 




3oj^my^^| 
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Fig. 10. Observed data at various tip-sample distances z. 

The PSF basis vectors were obtained from a standard PCA and the number of principal components (PC) in the 
PSF perturbation was selected as 4 based on detecting the knee in a scree plot. The same interpolation method as 
used in 1121 was adopted to account for unequal spatial sampling rates in the supplied data for the PSF domain 
and the image domain. 

In the PSF and image domains, along z axis, the grid in PSF signal space is 3 times finer than the observation 
sampling density, because the PSF sampling rate along z-axis is 3 times higher than the data sampling rate is. To 
interpolate this lower sampled data, we implemented a version of the Bayes MC reconstruction that compensates 
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for unequal projection sampling in z directions using the interpolation procedure of Dobigeon et al. |12| 
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(a) Ground truth synthetic virus (b) Semi-blind reconstruction of 
image obtained from data in De- the synthetic virus data. Only the 
gen et al (6). z-planes that have non-zero image 

intensity are shown.) 



Fig. 11. Results of applying the proposed semi-blind sparse image reconstruction algorithm to synthetic 3D MRFM virus image. 






(a) True PSF. 
Fig. 12. PSF estimation result. 



(b) Initial, mismatched PSF. 



(c) Estimated PSF. 



To demonstrate that the proposed 3D MCMC semi-blind reconstruction algorithm is capable of reconstruction 
in the presence of significant MRFM PSF mismatch, we first applied it to a simulated version of the experimental 



data shown in Fig. 10 We used the scanning electron microscope (SEM) virus image reported in Degen et al. |6| to 



create a synthetic 3D MRFM virus image, one slice of which is shown in Fig. 11(a) This 3D image was then passed 



through the MRFM forward model, shown in Fig. |12(a)| and lOdB Gaussian noise was added. The mismatched PSF 

1 2 = 0.0295. Fig. 



(Fig. 12(b) I was used to initialize our proposed semi-blind reconstruction algorithm. After 40 iterations the algorithm 

- tAuP from 0.7611 to Un- 



reduced the initial normalized PSF error 



I|k-oII 



11(b) 



and 



Fig. 12(c) show the estimated image and the estimated PSF, respectively. 

We next applied the proposed semi-blind reconstruction algorithm to the actual experimental data shown in Fig. 10 
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(a) MAP estimate in 3D and tlie estimated image on 6th plane, (b) Estimated (left) and nominal (right) PSFs. 11 tt^ — w^^ II = 

I'^ll ll^oli 

showing a virus particle. 0.0212. The difference between two is small. (Hard thresholding 

with level = raax{PSF) X 10^^ is applied for visualization. ) 

Fig. 13. Semi-blind MC Bayes method results and PSF coefficient curves. Az = A.3nm, pixel spacing is 8.3nm, X 16.6nm in x X y, 
respectively. The size of {x, y) plane is 498nm, X 531. 2nm. Smoothing is applied for visualization. 





(a) MMSE solution. Gray level (b) The pixel-wise square root of 
image intensity values range from the image variance. White color 
(black) to 7.34 X 10"^^ (white), indicates a high variance. Gray 

level image intensity values range 
from (black) to 3.29 X IQ-^^ 
(white). 

Fig. 14. The posterior mean and variance at the 6th plane of the estimated image (Fig. |13(a)|. 



for which there is neither ground truth on the MRFM image or on the MRFM PSF. The image reconstruction 
results are shown in Fig. [13] The small difference between the nominal PSF and the estimated PSF indicates that 
the estimated PSF is close to the assumed PSF. We empirically validated this small difference by verifying that 
multiple runs of the Gibbs sampler gave low variance PSF residual coefficients. We conclude from this finding 
that the PSF model of Degen et al. [6| is in fact nearly Bayes optimal for this experimental data. The blind image 



reconstruction shown in Fig. 13 is similar to the image reconstruction in Degen et al. [6] obtained from applying 
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the Landweber reconstruction algorithm with the nominal PSF. 

Using the MCMC generated posterior distribution obtained from the experimental MRFM data, we generated 

confidence intervals, posterior mean and posterior variance of the pixel intensities of the unknown virus image. 



The posterior mean and variance are presented in Fig. 14 for selected pixels. In addition, to demonstrate the match 



between the estimated region occupied by the virus particle and the actual region we evaluated Bayesian p-values 
for object regions. The Bayesian p-value for a specific region Ri having non-zero intensity is pv{Ri) = P({/fe = 
l}fee-Ri|y) where P is a probability measure and Ik is an indicator function at the kth voxel. Assuming voxel-wise 
independence, the p-values are easily computed from the posterior distribution and provide a level of a posteriori 
confidence in the statistical significance of the reconstruction. We found that over 95% of the Bayesian p-values 
were greater than 0.5 for the non-zero regions of the reconstruction. 

D. Discussion 

Joint identifiability is a common issue underlying all blind deconvolution methods, (e.g., scale ambiguity.) Even 
though the unicity of our solution is not proven, given the conditions that 1) span(K) = Kq + span(^ Ki) does not 
cover a kernel of a delta function, k ~ S{-), and 2) the PSF solution is restricted to this linear space of Kq, Ki's, 
the equation ( [TT] ) for the MAP criteria promises a reasonable solution that is not trivial such as x = y. Due to this 
restriction and the sparse nature of the image to be estimated, we can reasonably expect that the solution provided 
by the algorithm is close to the true PSF. A study of unicity of the solution would be worthwhile but is beyond the 
scope of this paper as it would require study of the complicated and implicit fixed points of the proposed Bayes 
objective function. 

Note that proposed sparse image reconstruction algorithm can be extended to exploit sparsity in other domains, 
such as the wavelet domain. In this case, if we define W to be the transformation matrix on x, the proposed 
semi-blind approach can be applied to reconstruct the transformed signal z = Wx. However, instead of assigning 
the single sided exponential distribution as prior for z, a double sided Laplacian distribution might be used to cover 
the negative values of the pixels. The estimation procedure for PSF coefficients, noise level, and hyperparameters 
would be the same. For image estimation, the vector h^ used in ( |22l ) would be replaced with the ith column of 
HW-i. 

VI. Conclusion 

We have proposed an extension of the method of the non-blind Bayes reconstruction in Dobigeon et al. fV2\ that 
simultaneously estimates a partially known PSF and the unknown but sparse image. The method uses a prior model 
on the PSF that reflects a nominal PSF and uncertainty about the nominal PSF. In our algorithm the values of 
the parameters of the convolution kernel were estimated by a Metropolis-within-Gibbs algorithm, with an adaptive 
mechanism for tuning random-walk step size for fast convergence. Our approach can be used to empirically evaluate 
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the accuracy of assumed nominal PSF models in the presence of model uncertainty. In our sparse reconstruction 
simulations, we demonstrated that the semi-blind Bayesian algorithm has improved performance as compared to the 
AM reconstruction and other blind deconvolution algorithms and non-blind Bayes method under several criteria. 

Possible extensions of the proposed method may include enforcing sparsity constraints on the result PSF and the 
eigenfunctions, by using sparse PCA type algorithms. Also, even with the selective sampling strategy that speeds up 
the sampling, the MCMC methods are slower than nonstochastic methods. This will be addressed in future work. 
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Appendix A 
Fast recursive sampling strategy 

In iterative MRFM algorithms such as AM and the proposed Bayesian method, repeated evaluations of the 
transformation T {k (A) , x) can be computationally difficult. For example, at each iteration of the proposed Bayesian 
myopic deconvolution algorithm, one must generate Xi from its conditional distribution / (x^jw, a, A^cr^,x_i,y), 
which requires the calculation of T {K,ii) where Xi is the vector x whose ith element has been replaced by 
0. Moreover, sampling according to the conditional posterior distributions of a^ and A^ (|26| and ( |23l ) requires 
computations of T (k, x). 

By exploiting the bilinearity of the transformation T {■,■), we can reduce the complexity of the algorithm. We 
describe below a strategy, similar to those presented in lfT2l App. B], which only requires computation of T(-, •) 
at most M x [K + 1) times. First, let Im denote the M x M identity matrix and u^ its ith column. In a first step 
of the analysis, the M vectors h,^ (i — 1, . . . , Af ) 

hf)=r(Ko,u,), (27) 

and KA4 vectors hf ^ (i = 1, . . . , M, k = l,...,K) 

hf)=r(vfc,u,), (28) 

are computed. Then one can compute the quantity T {k, x^) and T {k, x) at any stage of the Gibbs sampler without 
evaluating T (•, •), based on the following decomposition 

M K M 

T(/^,x) = ^x.hf)+^Afc^:r,hf). (29) 

The resulting procedure to update the ith coordinate of the vector x is described in Algorithm [3] below. 
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Algorithm 3 Efficient simulation according to / (x \w, a, ct , y ) 

At iteration t of the Gibbs sampler, for i = 1, . . . , M, update the ith coordinate of the vector 

T 



x(M-i) 



via the following steps: 
1: compute hi = hf' + Ylk=i '^kh.^i , 

2: 
3: 
4: 
5: 
6: 
7: 



(t) (t) (t-1) (t-1) (t-1) 



set T (/^, i|*""'^) = T (/^, x(*^^-i)) - xf -''h„ 
set e, =x-T('/«,if''"^^V 
compute fii, rjf and ly^ as defined in [6], 
draw x\ according to ( |2T| , 

set v(*^») - T-r^*) ^(*) ^(*) ^(*-l) ^(*-l)l^ 

set T (k, x(*^*)) - T f /«, i'^'''-^A + x(*'h,. 



Note that in step 7. of the algorithm above, T (k, x) is recursively computed. Once all the coordinates have been 
updated, the current T {k, x) can be directly used to sample according to the posterior distribution of the noise 
variance in (|26|. Moreover, this quantity can be used to sample according to the conditional posterior distribution of 
Xk in ( |23] l. More precisely, evaluating T {k (A^,) , x) in the acceptance probability ( p5] l can be recursively evaluated 
as follows 

T {k {XD , x) = T (/. (Ai*)) , x) - (aI*) - A*,) J2 ^»hf ^ • (30) 

i=l 

Appendix B 
Sparsity enforcing selective sampling 

Since we have estimated the 'overall sparsity', 1 — w, of x from ([T9|, we can expedite the sampling procedure 
of X by selectively sampling only significant portions of entire pixels of x. As a result, we expect (1 — w) x 100% 
of pixel domain of x to be zero, which will not need to be re-sampled. 

At time t, in order to approximate the quantile (1 — w) of {w,- }i=i,....M in pT| ) we evaluate the (1 — w) quantile 
value of the previously obtained {w^ }i=i....,M- This approximation accelerates the computation because the 
exact calculation of {wl }i=i m requires sampling of all x/s. Let q ^ quantile{{w\ ^ }i=i m, I — w) and 
Wthr — max((7, 1 — w). When w\ for x^ from ( |2T| ) is less than Wthr, then x[ is not updated or is set to zero. 
Because MCMC sampling is computationally expensive, especially for large size images, this suggestion can be 
restricted to the burn-in period to save computations. 

In our experiment, the selective sampling of x applied after 3 or 4th iterations produce equally good results 
compared to the conventional MCMC sampling methods, while reducing computation time by 30 — 50% for non- 
blind sparse reconstruction with a fixed PSF and by 10 — 30% for the semi-blind sparse reconstruction. 
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Appendix C 
Adaptive tuning of an acceptance rate in the random-walk sampling 

For an efficient sampling of Xk,k = 1,. . . ,K, from the desired distribution 7r(AA;) = / (A/c|A_fc,x, cr^,y), we 
need to properly tune the acceptance rate of the samples from the proposal distribution. A careful selection of a 
step size is critical for convergence of the method. For example, if the step size is too large, most of the iterations 
will be rejected and the sampling algorithm will be inefficient. On the other hand, if the step size is too small, 
most of the random walk moves are accepted but these moves are slow to cover the probable space of the desired 
distribution, and the method is once again inefficient. 

The transition density of Metropolis-Hastings sampUng is q{X^^\ A*^*-')acc(A^*-', A*^*^), where q{\'^*\ A**^*^) is the 
proposal density from A^*) and acc(A^*'', A**^*^) is the acceptance probability for the move from A'^*^ to A'^*^'^ Here 
we denote A/j by A without a subscript for simplicity. We set q{X^^\ A*^*^) to be a Gaussian density function of A'^*^*-', 
denoted by q{\'--*\ A**^*^) — g(A*'*-' — A*^*') = ^(A'^''*'; X^*\ s^) with a mean A^*^ and a variance s^, which produces 

a random walk sample path. Since q{-, •) is symmetrical, acCs(A(*\ A"^*^*^) = min ( 1, — , — , ) = 

\ 7r(A(*-')g(A(*-', A"^^*^) / 

min I 1, — , . j = P\(t)^\*(t), as derived in ( |24] i. Then the acceptance probability from a parameter value A*^*' 
is accs(A'^*') — J^,(t) (/(A*^*', A*(*^)acCs(A'^*\ A*'*^)dA*(*\ The acceptance rate with a scale parameter s, acting as 
a step size, can be expressed as accs = J^n{X)acCs{X)dX. 

We evaluate these integrations by using Monte Carlo methods; accg ~ — X)"=i acCs(A'-*') with A*-*' ~ 7r(A^*''), 
and acCs(A*^*-') w — ^"^-^ acCs(A'^*\ A*^*)) with A*^*) ^ q(A'-*', A**^*^). In practice, this empirical version of the 

integration value is evaluated as 

w 

ace, «— ^acc,(A(*\A*(*)), (31) 

after the burn-in period. Therefore we can evaluate the acceptance rate with s by averaging the Boolean variables 
of acc5(A'*\ A*'^*-'), t = 1. . . . , W, over a time-frame window of length W with realizations {A'*\ A*^*)}^. In short, 
we iteratively update s to achieve an appropriate acceptance rate, acCs, as described in Algorithm HI 

Algorithm 4 Tuning s in the Gaussian proposal density q{-,-) 

Select upper and lower limits acc^ and acc]^. At each time t — W, 2W, 3W, . . ., tune s via the following 

steps: 

1: Evaluate accg using ( |3T| l for the given time-frame window, 

s X c, if acCg > accH, 
2: Update s ^ { s ^ c, if acCs < accL, 
5, Otherwise. 



In practice, we fix the variance of the instrumental distribution at the end of a burn-in period. Consequently, the 
transition kernel will be fixed and this guarantees both ergodicity and stationary distribution. In our experiment, we 
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set a conservative acceptance range, ucch — 0.6, accL = 0.4, referring to ||29l , and W — 20,c = 4. This strategy 
can also be applied to the direct parameter estimation described in Appendix [Pj 

Appendix D 
Direct sampling of PSF parameter values 

As described in Section IT] in the MRFM experiments, the direct estimation of PSF parameters is difficult because 
of the nonlinearity of Kgen and the slow evaluation of Kg^niC') given a candidate value ^'. However, if Kg^n is 
simple and evaluated quickly, then a direct sampling of parameter values can be performed. To apply this sampling 
method, instead of calculating a linearized convolution kernel, k (A), we evaluate the exact model PSF, Kg^n (C)> 
in (|23| and ( p5] l. Also the proposed parameter vector C* correspondingly replaces a coefficient vector A* and the 
updated PSF is used in the estimation of other variables. This strategy turns out to be similar to the approach 
adopted by Orieux et al. in ll26l . 
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